clear	
set more off
if c(username)=="Scognamillo"	{
gl drive		"C:\Users\Scognamillo\Dropbox\FAO\Epic\Malawi\SPIA adoption\Antonio Scognamillo/IHPS"
gl tmp_panel	"$drive/tmp/josh"								
gl slm_do 		"$drive/panel_creation"
gl maps			"C:\Users\Scognamillo\Dropbox\FAO\EPIC\Spatial\Shapefile&coordinates\Malawi\MAP-Shape-Files"
gl graphs 		"$drive/graphs"
gl out			"$drive/out"
gl tables		"$drive/tables/raw"
								}
								
else if c(username)=="Antonio"	{
gl drive		"C:\Users\Antonio\Dropbox\FAO\EPIC\malawi\SPIA adoption\Antonio Scognamillo/IHPS"
gl tmp_panel	"$drive/tmp/josh"								
gl slm_do 		"$drive/panel_creation"
gl maps			"C:\Users\Antonio\Dropbox\FAO\EPIC\Spatial\Shapefile&coordinates\Malawi\MAP-Shape-Files"
gl graphs 		"$drive/graphs"
gl out			"$drive/out"
gl tables		"$drive/tables/raw"

}					 			






********************************************************************************
********************************************************************************
*********************************Maps*******************************************
********************************************************************************
********************************************************************************

set more off




use "$out\final_dataset.dta", clear
tempfile ssa_lab
lab save ssa_lab using "`ssa_lab'.do", replace


merge m:m district using "$tmp_panel/MW_SLM_district_label.dta",  nogen   keep(3)
replace district_label="Blantyre" 	if district_label=="Blanytyre"
replace district_label="Nkhata Bay" if district_label=="Nkhatabay"
replace district_label="Nkhotakota" if district_label=="Nkhota kota"
tempfile final_labeled
save `final_labeled'
*note Likoma and Neno district are missing in the using file

u "$maps/GADM/phdbADM1.dta", clear
rename NAME_1 district_label
rename id _ID
merge m:m district_label using `final_labeled',  nogen   keep(1  3)







preserve
drop if sq2==7
collapse (median)ssa_aez09 covrain_lt  mean_totrain_lt neg_shock_prob pos_shock_prob (max)sq2,   by(_ID district  time)


replace covrain_lt          		=round(covrain_lt	, 0.01)
replace mean_totrain_lt          =round(mean_totrain_lt	, 0.01)

do "`ssa_lab'.do"
lab val ssa_aez09 ssa_lab


spmap covrain_lt   using "$maps/GADM/phxyADM1.dta" if (time==1 | time==.), id(_ID) polygon(data($maps/DIVA-GIS/phxyWATER.dta)    fcolor(eltblue) oc(gs4)) /*line(data($maps/MWI_roads_xy.dta)   color(white ) size(medium ..))*/ fcolor(Oranges)oc(gs4) osize(thin) ndfcolor(none) ndocolor(gs4) ndlabel("no observations")   clnumber(5) legend(symy(*1) symx(*1) size(*1)) legorder(lohi) title(CoV rainfall) legend(pos(7)) saving("$graphs\cov_district", replace)
graph export "$graphs/cov_district.png",  as(png) width(800)   replace

spmap mean_totrain_lt   using "$maps/GADM/phxyADM1.dta" if (time==1 | time==.), id(_ID) polygon(data($maps/DIVA-GIS/phxyWATER.dta)    fcolor(eltblue) oc(gs4)) /*line(data($maps/MWI_roads_xy.dta)   color(white ) size(medium ..))*/ fcolor(Reds)oc(gs4) osize(thin) ndfcolor(none) ndocolor(gs4) ndlabel("no observations")   clnumber(5) legend(symy(*1) symx(*1) size(*1)) legorder(lohi) title(Historical mean rainfall) legend(pos(7)) saving("$graphs\mean_district", replace)
graph export "$graphs/mean_district.png",  as(png) width(800)   replace


spmap ssa_aez09   using "$maps/GADM/phxyADM1.dta" if (time==1 | time==.),  id(_ID) polygon(data($maps/DIVA-GIS/phxyWATER.dta)    fcolor(eltblue) oc(gs4)) /*line(data($maps/MWI_roads_xy.dta)   color(white ) size(medium ..))*/ fcolor(Reds)oc(gs4) osize(thin) ndfcolor(none) ndocolor(gs4) ndlabel("no observations")  clmethod(unique)   clbreaks(312 313 322 323) legend(symy(*1) symx(*1) size(*1)) legorder(lohi) title(Agroecological zone) legend(pos(7)) saving("$graphs\AEZ_district", replace)
graph export "$graphs/AEZ_district.png",  as(png) width(800)   replace

spmap neg_shock_prob   using "$maps/GADM/phxyADM1.dta" if (time==1 | time==.),  id(_ID) polygon(data($maps/DIVA-GIS/phxyWATER.dta)    fcolor(eltblue) oc(gs4)) /*line(data($maps/MWI_roads_xy.dta)   color(white ) size(medium ..))*/ fcolor(Heat)oc(gs4) osize(thin) ndfcolor(none) ndocolor(gs4) ndlabel("no observations")  clmethod(quantile)  clnumber(10) legend(symy(*1) symx(*1) size(*1)) legorder(lohi) title(Drought Risk) legend(pos(7)) saving("$graphs\drought_district", replace)
*graph export "$graphs/drought_district.png",  as(png) width(800)   replace
spmap pos_shock_prob   using "$maps/GADM/phxyADM1.dta" if (time==1 | time==.),  id(_ID) polygon(data($maps/DIVA-GIS/phxyWATER.dta)    fcolor(eltblue) oc(gs4)) /*line(data($maps/MWI_roads_xy.dta)   color(white ) size(medium ..))*/ fcolor(Blues2)oc(gs4) osize(thin) ndfcolor(none) ndocolor(gs4) ndlabel("no observations")  clmethod(quantile)  clnumber(10) legend(symy(*1) symx(*1) size(*1)) legorder(lohi) title(Flood risk) legend(pos(7)) saving("$graphs\flood_district", replace)
*graph export "$graphs/flood_district.png",  as(png) width(800)   replace

graph combine "$graphs\drought_district.gph" "$graphs\flood_district.gph" , /*title("`lab'" , size(medium)) subtitle(Rainfall Long-Run Mean and CoV, size(medium))*/  graphregion(color(white))
graph export "$graphs/drought_flood_district.png",  as(png) width(800)   replace

spmap sq2   using "$maps/GADM/phxyADM1.dta" if (time==1 | time==.),  id(_ID) polygon(data($maps/DIVA-GIS/phxyWATER.dta)    fcolor(eltblue) oc(gs4)) /*line(data($maps/MWI_roads_xy.dta)   color(white ) size(medium ..))*/ fcolor(Heat)oc(gs4) osize(thin) ndfcolor(none) ndocolor(gs4) ndlabel("no observations")  clmethod(unique)   legend(symy(*1) symx(*1) size(*1)) legorder(lohi) title(Soil Nutrients Availability Constraints) legend(pos(7) label(2 "No or Slight Constraint") label(3 "Moderate Constraint")label(4 "Severe Constraint")) saving("$graphs\nutrient_district", replace)
graph export "$graphs/nutrient_district.png",  as(png) width(800)   replace




restore






gl vars_panel   intleg_1   
foreach var of global vars_panel  {
recode `var' (.=0)
local lab `: var label `var''
preserve
collapse (mean) `var', by(ID_1 time) 
recode `var' (0=.)
gen `var'_perc = `var'*100
xtile q`var'_perc = `var'_perc, n(4)
tabstat `var'_perc, by(q`var'_perc) save
return list
mat qm= (r(Stat1) , r(Stat2) , r(Stat3) , r(Stat4))
sum `var'_perc
return list
mat minmax= (r(min) , r(max))
scalar min= minmax[1,1]
scalar q1= qm[1,1]
scalar q2= qm[1,2]
scalar q3= qm[1,3]
scalar q4= qm[1,4]
scalar max= minmax[1,2]
local min= round(0 )
local q1= round(q1 , 0.1)
local q2= round(q2 , 0.1)
local q3= round(q3 , 0.1)
local q4= round(q4 , 0.1)
local max= round(max, 0.1)
* play with round decimals if some district does not appear
cd "$drive"
local col1 "Reds" 
local col2 "Purples" 
local col3 "Greens"
local col4 "Purples"
local col5 "Greens" 
local col6 "Reds" 
local col7 "Reds"
local col8 "Oranges"
local i= `i'+1
spmap `var'_perc using "$maps/GADM/phxyADM1.dta" if (time==0 | time==.),  id(ID_1)  polygon(data($maps/DIVA-GIS/phxyWATER.dta)    fcolor(eltblue) oc(gs4)) fcolor(`col`i'') oc(gs4) osize(thin) ndfcolor(white) ndocolor(gs4)   clmethod(custom) clbreaks(`min' `q1' `q2' `q3' `q4' `max') ndlabel("no adoption or missing")   /*legtitle("% of plots")*/ legstyle(2) legjunction(-) /*ti("`lab'", size(medium))*/  legend(symy(*1.5) symx(*1.5) size(small)) legorder(lohi) subtitle("IHS 2010" , size(medium)) graphregion(color(white)) 
graph save "$graphs/`var'0",replace
spmap `var'_perc using "$maps/GADM/phxyADM1.dta" if (time==1 | time==.),  id(ID_1)  polygon(data($maps/DIVA-GIS/phxyWATER.dta)    fcolor(eltblue) oc(gs4)) fcolor(`col`i'') oc(gs4) osize(thin) ndfcolor(white) ndocolor(gs4)   clmethod(custom) clbreaks(`min' `q1' `q2' `q3' `q4' `max') ndlabel("no adoption or missing")   /*legtitle("% of plots") legstyle(2) legjunction(-) ti("`lab'" , size(medium))*/ leg(off) subtitle("IHPS 2013" , size(medium))  graphregion(color(white))
graph save "$graphs/`var'1", replace
graph combine "$graphs/`var'0.gph" "$graphs/`var'1.gph" , /*title("`lab'" , size(medium))*/ subtitle("`lab'", size(medium))  graphregion(color(white))
graph export "$graphs/`var'.png",  as(png) width(800)   replace
restore
}


global vars_ihps   mintill /*crop_rot*/ intleg_1 /*intleg_2*/ covercrop_any mulch  ///
		/*tree  intercrop_other  swc_mech swc_bio swc_any fert_org_use*/ 	///
		/*burn*/  ///
	

foreach var of global vars_ihps  {
recode `var' (.=0)
local lab `: var label `var''
preserve
collapse (mean) `var', by(ID_1 time) 
recode `var' (0=.)
gen `var'_perc = `var'*100
xtile q`var'_perc = `var'_perc, n(3)
tabstat `var'_perc, by(q`var'_perc) save
return list
mat qm= (r(Stat1) , r(Stat2) , r(Stat3) , r(Stat4))
sum `var'_perc
return list
mat minmax= (r(min) , r(max))
scalar min= minmax[1,1]
scalar q1= qm[1,1]
scalar q2= qm[1,2]
scalar q3= qm[1,3]
*scalar q4= qm[1,4]
scalar max= minmax[1,2]
local min= 0.0001
local q1= round(q1 , 0.0001)
local q2= round(q2 , 0.1)
local q3= round(q3 , 0.1)
*local q4= round(q4 , 0.01)
local max= round(max, 0.1)

cd "$drive"
/*
local col1 "Reds" 
local col2 "Greens" 
local col3 "Oranges"
local col4 "Purples"
local col5 "Reds" 
local col6 "Greens" 
local col7 "Oranges"
local col8 "Purples"
*/
local i= `i'+1
spmap `var'_perc using "$maps/GADM/phxyADM1.dta" if (time==1 |time==.),  id(ID_1)  polygon(data($maps/DIVA-GIS/phxyWATER.dta)    fcolor(eltblue) oc(gs4)) fcolor(Reds) /*fcolor(`col`i'')*/ oc(gs4) osize(thin) ndfcolor(white) ndocolor(gs4)   clmethod(custom) clbreaks(`min' `q1' `q2' `q3'  `max') ndlabel("no adoption or missing")   /*legtitle("% of plots")*/ legstyle(2) legjunction(-) ti("`lab'" , size(medium)) legend(symy(*1.5) symx(*1.5) size(small))  legorder(lohi) subtitle("IHPS 2013" , size(medium))  graphregion(color(white))
graph save "$graphs/`var'1", replace
graph export "$graphs/`var'1.png",  as(png) width(800)   replace
restore
}












